AMBIPOLAR DIFFUSION IN MOLECULAR CLOUD CORES 
AND THE GRAVOMAGNETO CATASTROPHE 



Fred C. Adams^'^ and Frank H. Shu^ 
ABSTRACT 

This paper re-examines the problem of ambipolar diffusion as a mechanism for the 
production and runaway evolution of centrally condensed molecular cloud cores, a pro- 
cess that has been termed the gravomagneto catastrophe. Our calculation applies in 
the geometric limit of a highly flattened core and allows for a semi-analytic treatment of 
the full problem, although physical fixes are required to resolve a poor representation of 
the central region. A noteworthy feature of the overall formulation is that the solutions 
for the ambipolar diffusion portion of the evolution for negative times {t < 0) match 
smoothly onto the collapse solutions for positive times {t > 0). The treatment shows 
that the resulting cores display non-zero, but sub-magnetosonic, inward velocities at the 
end of the diffusion epoch, in agreement with current observations. Another important 
result is the derivation of an analytic relationship between the dimensionless mass to 
flux ratio Aq = f^^ of the central regions produced by runaway core condensation and 
the dimensionless measure of the rate of ambipolar diffusion e. In conjunction with pre- 
vious work showing that ambipolar diffusion takes place more quickly in the presence 
of turbulent fluctuations, i.e., that the effective value of e can be enhanced by turbu- 
lence, the resultant theory provides a viable working hypothesis for the formation of 
isolated molecular-cloud cores and their subsequent collapse to form stars and planetary 
systems. 

Subject headings: MHD — methods: analytical — stars: formation 



1. INTRODUCTION 

Since molecular clouds are supported, at least in part, by magnetic fields, the removal of mag- 
netic fields represents an important component of the star formation process. In the most studied 
scenario, field removal occurs through the action of ambipolar diffusion, wherein magnetic fields 
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axe tied to the ionized component, which drifts relative to the more dominant neutral component 
of the gas (Mestel & Spitzer 1956, Mouschovias 1976, Nakano 1979. Shu 1983, Nakano 1984, Lizano 
& Shu 1989, Basu &: Mouschovias 1994). The neutral gas thus condenses inward toward a quasi- 
hydrostatic state, although perfect equilibrium is generally not reached. When the condensing gas 
becomes sufficiently centrally concentrated, the innermost regions of the structure begin to col- 
lapse dynamically onto a growing point-like protostar and eventually approach ballistic (free-fall) 
conditions. 

Because of its large dynamic range in space and time, the process is not easy to follow numeri- 
cally. The magnetic field and the forces that it exerts are vector quantities, so the relevant diffusion 
and dynamic equations are generally nonlinear, coupled partial differential equations in multiple 
spatial dimensions. The assumption of axial symmetry provides some simplification, but one still 
needs to deal with two spatial dimensions and time. In spite of these complications, the net result 
is physically simple: As magnetic fields diffuse outward, gas condenses inward to form a centrally 
concentrated structure that approaches pure power-law distributions of gas density, magnetic field, 
and other quantities (Nakano 1979, Lizano & Shu 1989, Basu & Mouschovias 1994). The goal of 
this paper is to demonstrate in a simple mathematical fashion how the asymptotic state is reached 
through nearly self-similar evolution toward a gravomagneto catastrophe, wherein an infinite cen- 
tral concentration is formally reached in finite time. After this point of catastrophe, set equal to the 
pivotal time t = 0, the system experiences true dynamical collapse to form a point-like star, which 
will be surrounded by a centrifugally supported disk due to the pre-collapse rotation of the cloud 
core. Note that this study does not include the eff'ects of rotation, so that the collapse solutions 
found herein represent the outer portion of the collapse flow; these solutions can then be matched 
onto inner solutions that include rotation and other effects (e.g., Cassen &; Moosman 1981, Terebey 
et al. 1984, Jijina & Adams 1996). 

Although molecular cloud cores experiencing ambipolar diffusion were identified as playing a 
dominant role in the formation of isolated low- mass stars two decades ago (Shu et al. 1987), recent 
observations indicate that ambipolar diffusion takes place more rapidly than the simple laminar 
description (e.g., Jijina et al. 1999). In addition, nonzero inflow velocities are often observed 
in starless cores (Lee et al. 2001, Harvey et al. 2002), which are assumed to be pre-collapse 
states. Contrary to popular perception, both of these properties can be accommodated within the 
picture of core formation via ambipolar diffusion. Indeed, one of the principal results of this paper 
is to calculate the nonzero, but sub-magnetosonic, inward velocity resulting from the ambipolar 
diffusion process. Despite this faster evolution, molecular cloud cores are well-deflned entities and 
not transient, turbulent phenomena (Lada et al. 2007). As a consequence of their linkage to strong 
magnetic fields, probably, these cores are also restrained from moving ballisticly through their 
parent clouds (Walsh et al. 2004). 

The rest of this paper is organized as follows. We present the formulation of the problem in 
terms of physical variables in §2, where we enforce axial symmetry and use a flattened approxi- 
mation. The following section (§3) outlines our approach to solving the resulting problem. We 
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first apply a similarity transformation, which converts the partial integro-differential equations into 
ordinary integro-differential equations. Since the diffusion time scale is comparable to or longer 
than the magnetosonic time scale needed to cross from the core center to the boundary where the 
core attaches to a common envelope, the solution to the induction equation itself requires a more 
complicated approach. In §4, we solve the zeoreth-order condensation problem for t < to describe 
the approach to the pivotal instant t = of gravomagneto catastrophe. As a simplification, we 
use the monopole (split monopolc) approximation for the gravitational (magnetic tension) forces, 
thereby transforming the integro-differential equations to ordinary differential equations that are 
solved by standard methods. In §5, we point out the shortcomings at large and small radii of the 
monopole approximation, and we generalize the approach by adopting various mathematical and 
physical fixes that show how the general physical problem possesses condensation solutions that 
connect smoothly to marginally critical envelopes. A central result of this section is the analytic 
derivation of a relationship between the dimensionless rate of ambipolar diffusion e and the dimen- 
sionless flux to mass ratio /o = Aq ^ of the central regions of the condensing core at the moment of 
gravomagneto catastrophe. In §6, we demonstrate how the runaway condensation that character- 
izes gravomagneto catastrophe transitions smoothly for t > to dynamically collapsing states that 
correspond to cores with accreting point-like protostars, i.e., the infall-collapse solutions that have 
been used widely in previous studies of star formation. In §7, we present a specific dimensional 
example to illustrate the typical astronomical characteristics of the entire process on both sides of 
the pivotal instant t = 0. We conclude in §8 with a summary of the astronomical implications of 
our results. Finally, in a series of Appendices A through F, we develop and extend various technical 
points encountered in the discussion of the text. 



2. FORMULATION 

The basic evolutionary equations for a flattened, self-gravitating, cloud core of surface density 
S and radial velocity u, threaded by a magnetic field with vertical component B^, are taken from 
the analysis of Shu & Li (1997) to be as follows. The equation of continuity is given by 

dE 1 d , , 

— + -—[ujuj:] = o. 1 

ot w ow 

The force equation is 

du du d „ 

where the acceleration produced by self-gravitation plus magnetic tension, -|- ^, is given by 
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and the kernel Kq is defined via 
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In equation (l2|), a is the gaseous isothermal sound speed, and Q provides the correction for the 
effects of the magnetic pressure (see Appendix A). FinaUy, the induction equation, which governs 
the evolution of the vertical component of the magnetic field threading the core in the presence of 
ambipolar diffusion, takes the form (see Appendix B) 



d 



S <! + —77— [tJ^uBz 

y at -co ow 



ZD dw 



27r-/C Si/2 



(5) 



where we have defined the radial component of the field at the upper vertical surface of the core by 

1 f°° 

Bt = — Ko{r/w)B,{r)rdr. (6) 
Jo 

The half-height zq appearing in equation ^ is defined by the assumed vertical hydrostatic equilib- 
rium (Appendix A). The quantities 7 and 1/C are, respectively, the usual drag coefficient between 
ions and neutrals and the height-averaged reciprocal coefficient for the ion mass-abundance (see 
Appendix B and Chapter 27 of Shu 1992). An attractive feature of the approach presented in 
this paper is that we can delay specifying the actual numerical value of the product 7C until it 
comes to specifying dimensional scalings appropriate to specific astronomical objects, as long as 
the combination of parameters given by equation (j20p below is a small number compared to unity. 



We define the dimensionless ratio A of mass per unit area to flux per unit area according to 

27rGi/2s 
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Appendix A derives expressions for @ and zq in terms of A for a magnetized singular isothermal 
disk, the form that our inner core approaches asymptotically at the moment of gravomagneto 
catastrophe. These relationships have the elegance of simplicity, and we adopt the approximation 
that the following expressions from Appendix A hold for all time, i.e.. 
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Combined with equations ([I]) , ([2]) , ([5]) , the relationships ([7]) and ([8]) give us a closed set of equations 
to solve for S, u, B^, A, 0, and zq. 



3. HOMOLOGY, SELF-SIMILARITY, AND ASYMPTOTICS 

In this section, we construct a similarity transformation to recast the problem in simpler form. 
First, we want to simplify the magnetic induction equation by using equation ([7]) and by making the 
(usual) approximation that the combination jC is a constant during the phase of molecular-cloud 
core formation. As a result, equation ([5]) takes the form 
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With this transformation, the definition of becomes 

L 



A(r) 



Similarly, the force terms now have the form 
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3.1. Basic Similarity Transformation 

With the form of the magnetic induction equation specified, we now transform from a {w, t) 
description to a (x, t) description, where we use the relations 

^(^'^) = 2'KG\t\ ^^^'^^ ' u{w,t) = av{x,t), B;,{zu,t) = -^^^^Pix,t) . (12) 

In the simplest type of transformation, self-similarity of the first kind (Barenblatt 1996), the func- 
tions a, V, and /3 introduced here would be functions of the similarity variable x only. In this case, 
however, we allow the functions to retain an additional time dependence to account for the fact 
that amibipolar diffusion occurs on a longer time scale than the runaway dynamics. Notice also 
that we have written the time variable in the coefficients with absolute value signs. We wish to 
mark the pivotal time t = as the moment of gravomagneto catastrophe, so that positive times 
correspond to the self-similar solutions of gravitational collapse onto a point-like protostar (Li & 
Shu 1997), whereas negative times correspond to the epoch of ambipolar diffusion in starless cores. 
The start of the ambipolar diffusion process thus corresponds to the limit t — > — oo, and the end of 
the ambipolar diffusion epoch corresponds to the limit t ^ 0~. 

With this choice of transformation, the dimensionless mass-to-flux ratio is given by 

X = Xix,t) = ^ (13) 
We also define its inverse, i.e., the dimensionless flux-to-mass ratio, 

/(x,t).^=|ii4. (14) 

X{x,t) a{x,t) 



The derivatives then take the forms 
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With this formulation, the equation of continuity is given by 
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The force equation then becomes 
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The induction equation can be written 
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so that e is a small dimensionless parameter of the problem (essentially the ratio of dynamical time 
to the diffusion time; see also Galli & Shu 1993, who denoted a similar inverse ratio as a large 
parameter %). In addition, the reduced radial magnetic field b is defined in terms of the integral 



1 /y\ 

b{x, t) = — ^0 ( - ) o-(y, t)f{y, t) ydy . 
X Jo ^x/ 



(21) 



If we use the traditional microscopic values of 7 = 3.5 x 10^^ cm^ g""^ s"^ and C = 2.0 x 10~^^ 

cni~3/2 gi/2 ^ggg Appendix B), we obtain e ~ 0.18. This small, but not very small, value of e 
allows for an illuminating, but not highly accurate, attack on the problem of molecular-cloud core 
formation and collapse. In practice, turbulence within the forming core may increase the effective 
diffusion coefficient by a factor of a several (Zweibel 2002, Fatuzzo &: Adams 2002, Heitsch et al. 
2004, Nakamura & Li 2005), which makes e a marginally small parameter. On the other hand, if 
cosmic ray fluxes are enhanced within molecular clouds (Fatuzzo et al. 2006), the value of e could 
be reduced by a factor of several. Thus, we anticipate that e might have large variations within 
molecular clouds, accounting in part for the wide range of observed core masses. The formal theory 
developed here allows a semi-analytical description of only those cores that form from regions where 
e ^ 1, but we anticipate that many of the physical insights gained from the formal analysis may 
carry over to the more general case even when e ~ 1. 



3.2. Solution by Iteration 

From many numerical simulations, we know that the effect of ambipolar diffusion in cloud cores 
is to try to redistribute the magnetic flux from the inner region, where the flux-to-mass ratio / has 
a relatively low, constant, value to the outer region where / has a relatively high, constant, value. 
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By relatively high, we mean typically / ~ 1; and by relatively low, we mean typically / ~ 1/2. 
Thus, / may vary only by a factor of 2 over a dynamic range in spatial scale of 10^, say, from 10~^ 
pc to 1 pc, whereas the volume density oc S/2zo over the same range of radii might differ by a 
factor of as much as 10^ (say, from 10^^ to 10^ molecules per cm^ - to be more precise, see Fig. [7]). 
As a consequence, it must be a good approximation to regard / to be a constant /o for calculations 
of the mechanical state of the most interesting parts of a condensing cloud core. 



To justify this conclusion mathematically, define a dimensionless measure of the time by 



(22) 



where to is an arbitrary unit of time used to make the argument of the logarithm dimensionless. 
To be definite, if we think of the core as having an outer boundary that connects to a common 
envelope at zuce, we may choose to to equal the time it takes a fast MHD wave traveling at speed 
to traverse the distance vjce- In numerical terms, this would typically make to ~ 10^ yr. 
With the definition (p2|) . equation ([19]) becomes 
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We proceed now to solve the governing set of equations by an iterative process. 



Begin by denoting solutions of the dynamical equations with / taken to be a fixed /o by the 
symbols ctq, vq, and b^. With the substitution of equation (I2ip . when b = bo, it is then trivial to 
show that the governing equation of continuity reads 



^ ^ ^"^"^^ ^ X ^ [xvo{x)do{x)] = , 



whereas the force equation becomes 
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(24) 

(25) 
(26) 



Once solutions for gq, vq, and 6o have been found, we can return to equation (p3|) and replace 
the factor /^/ \/l + /^ in the order e diffusion term by its zeroth-order approximation /o/\/l + /o • 
We can then obtain a better estimate for / by integrating the resulting linear partial differential 
equation for /, 
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where we have defined 

6o(x)^ r Ko(y)ao{y)^. (28) 

Note that the flux ratio f = fo has been removed from the definition of 60 and is now included 
in the leading coefficient. In principle, one could continue to iterate solutions of the flux-to-mass 
distribution from the induction equation with solutions of the surface density (from which we can 
get the magnetic field from the flux-to-mass distribution) and velocity field from the equations 
governing the mass and momentum flow of the fluid to obtain increasingly accurate numerical 
answers to the overall problem. In practice, we shall stop at the perturbative step (j27p . 



3.3. Homology and Self-similarity 

By introducing new scaled variables, we can transform the governing ordinary integro-differential 
equations for the zeroth-order dynamics into a universal form that is nominally independent of the 
numerical value of Aq = /q"^- Specifically, we adopt a scaling transformation of the form 

e = x/^/e'o, (29) 

v{0 = Mx)/^^^^, (30) 
a{O^Mx)il-f^)/^/e'o, (31) 

where the scaling coefficients are independent of x. Note that the flux ratio must obey the constraint 
/o < 1 (and hence A > 1) for the surface densities a and ao to be positive. The scaled forms of the 
equation of motion then become 

and 

5|._(i±^_^, (33) 

a at, ^ 

where the normalized force F is deflned by 

^(0 = -i r Koiv/O^iv) Vdv . (34) 



In the equations above, the discriminant D is given by 

V={^ + vf-l. (35) 



It is easy to see that the normalized equations (]32p and (133p are exactly what would have 
resulted if we had looked at the outset for self-similar solutions of the un-magnetized problem, 
/o = 0, Bq = 1. This mathematical homology explains the decades of confusion and controversy as 
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to what constitutes the proper "initial conditions" for the latter type of calculations (e.g., Larson 
1969, Penston 1969, Shu 1977, Whitworth & Summers 1985, Foster & Chevaher 1993, Andre et 
al. 2000). Our group (e.g., Shu, Adams, & Lizano 1987, Lizano & Shu 1989, Li & Shu 1996) has 
long maintained that the pivotal instant t = represents, not an "initial condition" where singular 
conditions are reached at the origin, but rather a transitional instant between an extended period, 
t < 0, of magnetic evolution through flux loss via ambipolar diffusion, and another period, t > 0, 
of dynamical collapse, infall, and star plus centrifugal disk formation. After some initial debate, 
the group led by Mouschovias has come to the same conclusion (see, e.g., Tassis &: Mouschovias 
2005). This paper then provides the mathematical justification for the latter point of view, and 
it supplies the means to select from the wealth of "extended-contraction/runaway-condensation" 
solutions for t < advocated first by Hunter (1977) as worthy alternatives to Shu's (1977) choice to 
start at t = with singular isothermal systems at rest, after what Shu argued would be a period of 
subsonic evolution to reach such a state. The corresponding static starting state here reads: v = 0, 
a = 1/^, and F = —1/^, which provide exact solutions of the equations ([32]) . ([33|) . and but 
not exactly those that we want here. 

The critical point of the flow occurs where V = 0. In order for the flow to pass smoothly 
through the critical point the right-hand sides of both equations ([32]) and ([33j) must vanish 
where D = 0. This requirement defines two conditions, which act to fix the values of v and a at 
the critical point ^ = Using L'Hopital's rule, we can integrate inwards and outwards from 
We require that the inward integration satisfies the inner boundary condition v = at ^ = 0. Note 
that only one value of ^* can satisfy this constraint. With the critical point thus specified, the 
outward integration from the same point = ^* produces the asymptotic behavior a A^~^ and 
V — > Voo as ^ — > cxD. 



3.4. The Flux to Mass Distribution and Intermediate Asymptotics 

In reduced and scaled variables, the equation (j27p can be written 



9 In r 

where we have defined 



e^2(o 



^ [C'^imo] , (36) 



, . (37) 

Note that the dependence of the overall scaled problem on the parameters e and /o enters explicitly 
only in equation (|36p through e. In what follows, we shall see that the proper formulation of an 
initial- value problem for the solution of equation ()36l) will connect e and /o, i.e., that the flux-to- 
mass ratio of the central-most regions of a condensing cloud core depends on the dimensionless rate 
of ambipolar diffusion as measured by the parameter e. 
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The Unear partial differential equation (j36p can be attacked by the method of characteristics: 
df d 1 

— = IhfiE) on the trajectory — (— Inr) = -— , (38) 

where 

is regarded as a known function of ^. Thus, we may define the formal integral, 



NiO = J^m)dC and ^0 = ^ 
and we write the solution to equation (j36p as 



(40) 



/(e,r) = /(l,ri) + eiV(0, (41) 

where ^ is the present position at the present time r connected to a past (or future) position 1 at 
the time ri on a Lagrangian trajectory that reads in similarity coordinates: 

r(0 +lnr = Inn. (42) 



With equation ([42]) giving ri as a function of ^ and r, equation ()4ip yields the general solution 
for the advection-diffusion equation (|36p where the term — eA^(^) gives the effect of the ambipolar 
diffusion relative to a comoving observer and /(l,ri) gives the effect of advection if we follow the 
fluid motion assuming field freezing. 

To illustrate the behavior of /(l,ri), we first note that the position ^ = w/Q^^'^atoT = 1 lies 
just outside the origin = as r — > 0, whereas the same = 1 position lies at a great radial 
distance w in the limit r ^ oo. Thus, we have the generic behavior: 

/(l,0) = /o, /(l,oo) = l, (43) 

if the supercritical core connects to a marginally critical common envelope. In §5.3, for reasonable 
core models, we shall find that N(oo) is small compared to unity. [See Fig. 4 and note N{oo) = 
No{oo) — A'^o(l)-] On the other hand, note that if f(^) were zero, T(^) would equal In^, and the 
characteristic trajectory would simply follow a line defined by In(r^) = const, i.e., a line of constant 
= Tu/QoatQ or constant zu (because fluid elements are not moving if v were zero). Although the 
inflow velocity, f (^), is not zero in our problem, it becomes a constant at large where the term 
^ + f (C) is dominated by ^. The relationship between ri and and r is then given by 

Ti w ^r. (44) 

Thus, in the limit r — > with ^ 1, when the moment of gravomagneto catastrophe is approached, 
the flux to mass distribution as given by equation (j4ip has the approximation /(^,t) ~ fi^,(,T) 
and assumes all intermediate values between /(1,0) = /o at small oc to /(I, oo) = 1 at large 
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oc w. In other words, during the runaway phase of core condensation, the function /(^, r) 
"freezes" with a profile that is a function only of the Lagrangian coordinate (which could be taken 
to be the enclosed cynlindrical mass) varying monotonically from /o at the core center to unity at 
the outer core boundary. 

Consider now a position ^ <C 1 at small but finite r > after runaway condensation is in 
progress (which occurs roughly at r ~ 1; see §5.3), with the spacetime point (^,t) being connected 
to an initial pair (1,ti) near the outer core boundary, i.e., where /(C,t) /q and /(l,ri) « 1. 
Equation (|^T|) then requires 

- eiV(O) = 1 - /o. (45) 

Together with the definition ()37p . equation ()45p provides us with an eigenvalue relationship between 
/o and e. In other words, for given e, runaway condensation occurs when ambipolar diffusion has 
produced a central flux to mass rato that satisfies 



(l-/o)v^rT2:f = e (46) 

where A^o(l) = -A^(O) = /q^ AA(0 di (see Fig. 4 in §5.3). 

Unfortunately, these properties of the general solution depend on the seemingly innocuous 
assumption that Mi^) is integrable at ^ = and cxd. However, as shown in the following sections, 
at small ^ the functions ct(^), f(C), and approach the forms 

^(O^^o, ^(e)^-e/2, F(C)^F'(0)e, as e^O. (47) 
Thus, for small ^, AA(^) behaves as 

M{i) = ^ for (48) 

In these circumstances, N(^,t) will diverge as [4:F'(0)/ao]ln(^ as ^ — > 0. The divergence arises 
because in the derivation for the average value of in equation (B3), we have set dB^^/dz 
equal to (-B+/2o)sech^(2:/2;o). This approximation is valid away from the origin, but at the origin, 
dB^/dz is doubly small, because not only B^ oc B^ itself is small, but the magnetic field is vertical 
near the origin so d/dz is also small. The replacement of dB^/dz by something proportional to 
B^/zq oc o'(^)F(^) accounts for the first effect, but not the second. 

As a related point, the current density oc {dB-^/dz—dBz/dw) is dominated at the origin, not by 
dB^/dz, but by dB^/dvj, implying that the Lorentz force there comes mostly from the gradient of 
the "magnetic pressure" —d{B'^/8Tr)/dw rather than from the "magnetic tension" {Bz/4:7r)dB.^/dz. 
This dominance is evident in that it is the pressure (gas plus magnetic) that decelerates the inflow 
to rest at the origin, not the tension (see eq. [73|) . The transition in roles of the tension versus 
the pressure when one moves from the disk of the core to its central regions has implications 
for the ambipolar diffusion that occurs near the origin, which the unadulterated diffusion term 
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M in equation (j39p does not treat correctly. Indeed, near the origin, the diffusion term involves 
the second derivative d'^B^/ dw"^ that translates into a term proportional to d'^{af)/dS,'^. In a 
rigorous discussion, we would examine the central regions anew and provide a proper matching of 
the solutions there with those applicable for the more highly flattened regions of the cloud core 
developed here. The extra radial derivative that appears via d'^{af)/d^'^, multiplied by a small 
coefficient e, makes possible asymptotic matching of the inner solution to an outer or intermediate 
solution. In particular, it would be possible to invoke an extra boundary condition, say, df /d(, = 0, 
at the origin to ensure that / is well-behaved at the origin. A formal attack along these lines would 
involve singular perturbation theory, coupled with the introduction of multiple length and/or time 
scales (Bender &: Orzag 1978). 

In other words, our problem is really one of intermediate asymptotics (Barenblatt 1996). A 
proper treatment would asymptotically match the intermediate core on both an outer scale to a 
common envelope, where the assumption of gravitational contraction breaks down. It would do the 
same on an inner scale to the central region, where the assumption of a flattened configuration is 
invalid. For the sake of physical clarity, we forego such a formal study in this paper, and join the 
intermediate core to a common envelope only in terms of its magnetic connection and not in its 
the mechanical considerations. We also incorporate the region near the origin into the intermediate 
analysis by flxing the problem presented above, as well as one that will appear later, with simple 
procedures that are physically rather than mathematically motivated. We defer to §5.3 therefore 
our prescription for making AA(^) regular at the origin. 



4. SHEET CONDENSATION SOLUTION 

4.1. Gravity Formulation in Terms of the Enclosed Mass 

In order to obtain an approximation to the condensation solution for t < 0, we can make 
the assumption that the potential is given by the monopole associated with the enclosed mass. 
Although it is possible to develop this approximation as the first term of a general multipole 
expansion (see Appendix E), we prefer a motivation based on physical intuition. We start by 
defining the cylindrically enclosed mass M{w, t) in dimensional units: 

M{w,t)= / 2TTrdrT.{r,t) . (49) 
Jo 

The statement that the enclosed mass is conserved if we follow the motion of mass annuli, 

dM dM ^ 
at ow 

may be regarded as an integrated form of the continuity equation ([1]). If we use the dimensionless 
variables defined by equations (fl^ . then equations (jlUj) and (fSOj) become 

M{'uj,t) = —\t\m{x) where fh{x) = / a xdx , (51) 
^ Jo 
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din 

-m+{x + v) = 0, (52) 

ax 

where we have assumed the case where t < 0. We employ the same scahng transformation as 
before (see eqs. [29] ~ [31]) to express the equations in terms of the variables ^, v, and a. The 
dimensionless enclosed mass now becomes 

e3/2 .5 
rh{x) = — ^—-^m{S^), where nT-iO = / cr ^d^. (53) 
1 - Jo ^0 

The differential version of the last equation, 

^ = ia, (54) 
may be combined with the scaled version of equation ([52]) . 

diTTl 

m={i + v)— = i{i + v)a, (55) 
to express the monopole approximation for the force in the form 

F{0 = -^ = -\{i + v)<y- (56) 

As a check, notice that ^ times the equation of continuity, written in the usual fashion as the 
scaled version of equation (f24l) . 



is simply the derivative of equation (j55p with respect to ^. With the force F reduced to a local 
expression, equations (j32p and (j33p become the following coupled set of first-order, nonlinear, 
ordinary differential equations: 

p|=(i±^ II + (58) 
P^ = (f + -({ + ")] , (59) 



where the discriminant T> is given by equation ()35p . 



4.2. Critical Points 

At the critical points P = 0, so the right hand sides of equations (f59l) and (f58l) must also 
vanish. This condition thus determines the value of the density field at the critical point, 

cT(e*) = l. (60) 
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If we expand around the critical point, the leading order corrections have the forms 

£, = £,* + SS, , V = + vi6S, , and cr = ct* + aiS^ . (61) 

Using these expressions in the equations of motion and keeping only the leading order terms, we 
can find the field derivatives at the critical point: 

v, = -^±±- [(^, - If + 1] and ai = ^ {e* - 2 T [{C* - if + l] '^'} ■ (62) 

Appendix D generalizes this procedure for arbitrary forms of F(^). 



4.3. Limiting Forms 

In the limit ^ — ^ oo, the force equation allows for an asymptotic solution for the surface density, 
namely 

a = - , (63) 

where A is a constant. With this form for the surface density, the asymptotic behavior for the 
velocity can be found, 

v{i) ^ ^^oo + , (64) 

where v^o is a constant. The second correction term goes to zero as ^ ^ oo, but the falloff is slow. 
As a result, the inflow velocity will in general be nonzero but small even at the pivotal instant of 
gravomagneto catastrophe (Allen et al. 2003, Fatuzzo et al. 2004); this result, in turn, implies 
that the subsequent infall rates will be larger than for cases where the starting states are in exact 
hydrostatic equilibrium. 

In the limit ^ ^ 0, we want to enforce the boundary conditions 

v — ^ and cr — > (To = constant . (65) 
The solution for the velocity field near ^ = has the dependence 

<K) = -\i- (66) 
With this limiting expression for u, the density field takes the form 

Notice that this particular function becomes cr = 2/^ at large values of the scaled similarity variable 
^, whereas the asymptotic limit of the density field has the form a = Aj^. 

Unfortunately, the behavior near the origin is less than perfect in a model of starless cloud cores 
as flat sheets. The monopole associated with a sheet produces a net scaled force (gravitational plus 
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magnetic tension) that approaches the form F — (To/2 = —2.524 as ^ ^ (see eq. [56] and below). 
The origin of this behavior rests with the unsealed reduced force behaving as -F = — (1 — fQ)fhQ/x^ 
with rfiQ = aox^/2 if the surface density approaches a constant o"o in the limit x ^ Q. In other 
words, with constant surface density, the mass enclosed within a cylindrical radius scales as the 
square of that radius, which cancels the inverse square law of Newtonian gravity for monopoles (and 
Amperian magnetism for split monopoles). Hence, unlike the classic freshman physics problem of 
the gravitational field inside a sphere of uniform volume density, the corresponding force in an 
axisymmetric flat sheet of uniform surface density does not go to zero, even when we approach the 
origin! 

However, by symmetry considerations, the radial force must vanish right at the origin for 
any axisymmetric mass and current distribution that is regular there. Appendix E shows that 
a jump in the physical behavior as one changes from being at the origin to being slightly off it 
is generic to any order in a general multipole expansion of an axisymmetric sheet. On the other 
hand, incompletely flattened, magnetized, isothermal, disks/toroids have scaled reduced half- height 
Co and dimensionless volume density of the form a = (cr/2Co)sech^((^/(^o) (see Appendix A). If the 
volume density is regular at ^ = (rather than diverging as which applies only at the pivotal 
instant t = 0), then Foc^^Oas^— >0 (see Appendix F). Accounting for the finite thickness 
of actual molecular-cloud cores thus cures the unphysical situation at the origin. As a result, a 
solution that enforces physical boundary conditions at the origin based on the fluid's reaction to 
a sheet monopole is clearly blemished. In the interest of obtaining practical and useful results, 
however, we defer further discussion of this imperfection until §5. 



4.4. Sheet Monopole Solution 

With the preliminaries in place, we can now find the critical points, and integrate both inward 
toward ^ = and outward to large ^ to find the solutions for the reduced and scaled density 
and velocity fields (in this monopole approximation). In the usual case, these two coupled ODEs 
would require two boundary conditions to specify a solution. In this setting, equations ([65]) supply 
the inner boundary conditions on v and a, but we do not know the correct value (Tq to enforce. 
However, this problem contains an additional constraint, namely, that the flow must pass smoothly 
through the critical point (specifically, the fluid fields v and a must be continuous at but 
their derivatives need not be). Since each starting value of do would lead to a different value of 
a at the critical point, only one value do allows for smooth flow. To find this value, we start the 
integration at a possible critical point and integrate inwards toward = 0. By requiring that the 
solution satisfy the inner boundary condition on v, we can iterate the starting point until we find 
the correct value of the critical point. With this value specified, we then integrate outward from the 
critical point. With no further quantities to specify, this integration thus determines both A and 
Voo ■ The resulting solution is shown in Figure [1] (along with solutions from the following section) . 

The solutions for v and a follow the limiting forms found analytically in §4.3. From the 
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Fig. 1. — Reduced fluid fields for the pure monople, monopole plus quadrupole, and full sheet 
solutions plotted, respectively, as the solid, dotted, and dashed curves. The solutions for reduced 
and scaled surface density a{^) and velocity v{^) are given, respectively, by the upper and lower 
panels, with the limiting values a ^ ao and v^Oas^^O; A/( and t; — > t^oo as ^ — > oo. 
The parameter values are tabulated in Table 1. 
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numerical solution, we can find the values of the parameters appearing in the analytic forms: 
= 1.294, (To = 5.03, A = 1.40, and — foo = 0.495. We can combine the inner and outer limiting 
forms for the density profile to construct an approximate solution, 

that agrees well with the numerical solution with an RMS error over the range of ^ shown in Figure 
[T]of only ~ 1%. Similarly, we can fit the velocity field with the form 

C + 2|Voo| 



where this result agrees with the numerical solution with an RMS error of 



70. 



One of the most important quantities resulting from this calculation is the value of the nonzero 
velocities in pre-collapse cores, i.e., cores that have not yet produced protostars at their centers. 
The physical value of this inward speed is given by Uoo = avo(^ ^ oo) = aVooVOo- Since the flux 
to mass ratio /o must lie in the range < /o < 1, the correction parameter Qq is confined to the 
range 1 < Go < 3/2. As a result, the head-start velocity is constrained to lie in the range 

0.495 < ^ < 0.606 , (70) 

a 

where a is the isothermal sound speed (and the minus sign denotes inward velocities). These values 
are consistent with, although perhaps slightly smaller after projection, than the extended infall 
velocities observed in starless molecular-cloud cores (e.g., Lee et al. 2001, Harvey et al. 2002). 

Another important physical quantity in this problem is the size of the region of nearly constant 
density in the core center. The similarity solution (Fig. [T|) shows that the surface density a has 
zero slope in the center and steepens to the form a = A/S, in the outer regime. We can thus define 
the outer boundary of the core region to be the location where the index p = —{i/<j){da/d^) = 
1/2. Using the approximate form given by equation (I68p . we find that the outer boundary of the 
core region occurs at ^1/2 ~ 0.313. The physical location of this boundary is given by 

tui/2 w 0.313 a |t| ^/Wq, (71) 

where the column density falls to 0.519 of its central value. 

To summarize, before t = 0, condensing, magnetized molecular cloud cores display a finite re- 
gion of nearly constant central surface density that makes them mimic static Bonnor-Ebert spheres. 
In actuality, however, the surrounding regions of the core are in a state of extended contraction at 
a significant fraction of the isothermal sound speed a. In the limit t — > 0~, this central region loses 
its finite extent and the core attains a pure-power law configuration, S oc ro^^. that corresponds 
ideally in three dimensions to a flattened singular isothermal toroid, p oc r'''^R{0) in spherical polar 
coordinates. The core subsequently goes into true dynamical collapse onto a central protostar. 
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5. 



BEYOND THE SHEET MONOPOLE APPROXIMATION 



Before we consider the collapse solution for t > (see §6), we consider the three shortcomings 
in our treatment of the gravity of a flattened core by the sheet monopole approximation. 

[1] For a flattened core, the monopole approximation represents only the first term in a more 
general multipole expansion (see Appendix E). These higher- multipole terms convey two types of 
corrections: Inner multipoles correct for the fact that matter in a flat sheet interior to the field 
point ^ is (on average) closer to ^ than if the enclosed mass were placed at the core center. Outer 
multipoles correct for the gravitational pull of the matter in the sheet outside of the field point 
^. The physical description for the action of currents is more complicated because we have used 
a scalar potential rather than the vector potential to describe the magnetic field (see Li & Shu 
1997), but the consequence for the magnetic tension is the same except the tension force acts in 
opposition to self-gravity. 

[2] The aspect ratio, which is given by 



cannot be small compared to unity in the limit ^ ^ 0. Indeed, even at large ^, zq/tu = 2(1 — 
/o)/A(1 + 2/q), which equals 1.4(1 — /q)/(1 + 2/q) for the sheet monopole solution, and is not small 
unless /o is close to unity. The assumption that the cloud core is highly fiattened is egregiously vio- 
lated in the central regions, precisely where the surface density profile flattens instead of continuing 



which has a value 1.22 for a typical /o = 1/2. 

[3] The regions at large ^ are not fully isopedic with constant f = fo (see §5.3). If / is an increasing 
function of ^, then the force of magnetic tension becomes increasingly strong relative to the force of 
self-gravity, instead of maintaining a constant ratio (with opposite signs), as is true in the inner core. 
Growing magnetic support against self-gravitation as the envelope is approached will presumably 
also reduce the induced inflow velocities. 

We now discuss how these shortcomings may affect a peculiar aspect of the sheet monopole 
solution obtained in §4. We found that the flow properties are completely deflned by the behavior 
of the gas near the origin ^ = 0. In particular. Appendix C proves that the velocity profile is 
monotonic, which implies that if inflow occurs at any point in the self-similar system, then (for 
t < 0) the flow must pass through a critical point where v = 1 — ^,, and approach the form v = —(,/2 
near the origin. The latter behavior during the cloud-core condensation stage has nothing to do with 
gravitational minus tension forces. It represents the deceleration of an initially inwardly directed 
velocity by the pressure forces. Independent of how the force F behaves, as long as it does not 
diverge at the origin, equation (j32|) implies at small ^, 




(72) 



inwards as o" = A/S,. For the monopole solution, at = ^1/2 we have zq/w = 2.44(1 — /q)/(1 + 2/q) 



, ,.dv 1 , ^ , 



(73) 
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where the —1 remaining in the discriminant comes from the pressure gradient. The above equation 
has the solution v = —£,/2 if the velocity v vanishes at the origin. 

In itself, the above result is not particularly ominous. But smooth passage through the critical 
point also determines the solution that is reached at asymptotic infinity, in particular, the values of 
the head-start velocity Voo and the surface-density coefficient A. How is it possible for the conditions 
near the center of a condensing cloud core to dictate how the core connects at asymptotic infinity 
to the cloud envelope? Isn't that putting the cart before the horse? In principle, if there is enough 
time and the inflow is sub-magnetosonic, as is the case with the monopole solution, then one 
could imagine the central regions to have a magnetohydrodynamic influence on the outer regions. 
However, we are not guaranteed such simple behavior in every circumstance, and naive treatments 
of the core gravity can produce super-magnetosonic condensation speeds (see below). Wouldn't 
such solutions be unstable to shock formation as the pressure forces attempt to bring the inflow to 
a halt at the core center? (Compare this question with previous criticism [Shu 1977] of the Larson- 
Penston solution, which is over-dense and supersonic by even larger margins, and the placement of 
the Larson- Penston solution in the context of champagne outflows [Shu et al. 2002].) 

Slow core condensation by ambipolar diffusion avoids the above paradox. As flux is lost from 
the central regions to the outer regions above (to reach typical values of /o ~ 1/2), the inner core 
begins a stage of extended contraction at a fraction of the magnetosonic speed ~ as is seen in both 
numerical simulations (e.g., Basu & Mouschovias 1994) and observations (e.g., Lee et al. 2001, 
Harvey et al. 2002). If the leakage of the flux is slow, as it must be because the envelope of a 
typical molecular cloud is too well-ionized to allow rapid motion of neutrals past ions, then the 
condensing core never loses so much support that its surface-density coefficient in the outer parts 
becomes substantially greater than the equilibrium value A = 1. Without large over-densities A—1, 
it is not possible to generate large head-start velocities Voo - unless one has an over-idealized force 
calculation for -F(^) in the central regions by assuming the region is highly fiattened when it is not. 

With the above comments in mind, we extend the fiattened monopole approximation by two 
different methods. In the first method, we retain the sheet approximation but compute the force 
in its full form (see Appendix E): 



which is mathematically identical to equation (|34p . The difference surface density A/r] — cr{r]) is 



integral may be truncated at a reasonable upper limit without compromising numerical accuracy. 
Appendix D gives a formal description of the solution procedure for such arbitrary forms of -F(^). 

In the second method, we retain the interior monopole approximation, but compute the force 
in a modified form. We begin by defining the reduced scaled half-thickness 




(74) 



everywhere positive but rapidly goes to zero much outside of the central core r] ^> ^1/2- Hence, the 




(75) 



- 20 - 



We then replace ^ in the denominator of equation (j56p by y^C"^ + CoiO the heuristic basis that 
the latter is a truer measure of the distance between the field point and a typical interior source 
point in an incompletely flattened cloud core. The softened-monopole force -F(^) now reads 

FiO = - (76) 

ve+cm 

Note that equation (j76p produces a force law that is proportional to ^ for small ^ where v ~ — 
and a ~ fio. Notice also that that the homology with unmagnetized systems has disappeared; cases 
with different /o produce different reduced scaled variables. In partial compensation, forces of the 
form ()76p are directly integrable by a slight modification of the method discussed in §4. Finally, 
notice that the implied "enclosed mass" no longer corresponds to the cylindrical value because in a 
quasi-spherical geometry, it is more appropriate to consider a central region with constant volume 
density, which still integrates to a constant surface density in the complete vertical direction, but 
does not behave at small ^ as a sheet of matter and current. 



5.1. Full-Gravity Corrections in a Sheet 

The case when -F(C) is the full gravity of a sheet is shown as the dashed line in Figured! Note 
that the solution has basically the same form as the sheet monopole solution (solid curves), but 
exhibits somewhat different values of the defining constants. The critical point in the flow shifts 
outward to ~ 1.714; the central density do becomes somewhat smaller 1.98, while the asymptotic 
density coefficient A increases to 3.43, resulting in the larger (supermagnetosonic) head-start speed 
—Voo ~ 1.67. For reference, the intermediate case for the monopole plus quadrupole corrections 
in a sheet geometry is plotted as dotted curves. We may attribute the full-gravity and monopole- 
plus-quadrupole results to net inward gravitational fields that have increased in the outer regions 
and decreased in the central regions relative to a pure monopole. 



5.2. Softened Monopole Gravity 

Unfortunately, the approach of the previous subsection has its own difficulty in the central 
regions because the approximation that the disk is geometrically thin leads to an unphysical, non- 
vanishing, force -F(^) near the origin before a protostar has formed there. Equation (j74p shows 
explicitly how tricky it is even to conclude that -F(^) equals a constant at the origin rather than 
diverging as 1/.^ when we model a cloud-core with a power-law envelope and a non-singular center 
as a flat sheet. 

The behavior at asymptotic infinity for the softened monopole is similar to the sheet monopole, 
but it is somewhat different near the origin. The velocity field has the same form as before, i.e., 

v^-^^ as e^O. (77) 



- 21 - 



The density field approaches the form 



(0 = ^ [l - (1 - co/ag) e-€'/4 



1/2 



(78) 



where we have defined 



co = (l-/o')/(l + 2/o'). 



(79) 



Note that, instead of a finite value, the gas pressure gradient is now zero at the origin, reminiscent 
of quasi-hydrostatic equilibrium. 

Figure [2] shows the different approximations for the force prescription used in this paper. The 
dashed curve shows the force law for the case of full sheet gravity. Since its force changes sign at 
small ^ (the net force as ^ — > is constant but directed outward)^ we cannot display the innermost 
region in a log- log plot. The dotted curve shows the force profile resulting from the sheet monopole 
approximation for the same surface density and velocity profile that resulted from the case for the 
full sheet gravity. This normalization is needed to make a fair comparison of the three techniques. 
Note that the sheet monopole force approaches a constant value (now directed inward) in the limit 
^ ^ 0, which leads to the difficulties discussed earlier. Finally, the solid curve shows the force 
profile calculated using the same rules but with the softened monopole approximation for the case 
/o = 0.5. By construction, its force law has a linear form at small ^, vanishing at the origin, and 
then joins onto the standard profile at larger values of ^. 

To study the possible range of head-start velocities, and to assess the meaning of the artificially 
inflated values resulting from the anomalies of all forms of sheet gravity at the origin, we use the 
modified monopole prescription of equation (j76p to find additional condensation solutions. The 
resulting solutions can also be characterized by the parameters fio, A, and — Woo; which are tabulated 
in Table 1 for varying values of /q. Figure [3] plots the associated functions C7(^) and i'(C)- Notice 
that as the dimensionless flux-to-mass ratio /o decreases, the degree of softening becomes larger, 
the effective scale height of the density distribution in the radial direction becomes larger, and the 
critical points of the flow move outward. The flow solutions in the outer regime are specifled by 
the parameters Voo and j4, which grow larger with decreasing /q: The softer gravity allows for more 
extended density profiles and hence larger A, and pulls inward less strongly to allow for greater 
head-start velocities Voo (see Table 1 and Figure [3]). The important point, however, is that the 
head-start speed —Voo in each case is suh-magnetosonic (i.e., \voo\ < !)• 

Notice also that in the limit /o ^ 1, the softening parameter ( ^ 0, and we recover alge- 
braically the unsoftened monopole approximation. As a result, the /o ^ 1 solution agrees with the 
original monopole solution (see Fig. [3] and Table 1). This case is to be regarded as the limiting case 
where /o — > 1 from below. In this limiting procedure, the weak effective gravity is compensated by 
the large dimensional surface density, and thus produces roughly the same collapse time scale (see 
eqns. [88] and [99]): M^ore/M = {A/mo){ujce/Q'^^'^a) for cores with different values of /q. 

Since a realistic treatment should include both of the effects discussed in this subsection and 
the previous one, we expect the actual results for uo. A, and —Voc to be given by a combination 
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Fig. 2. — Comparison of the force profiles calculated from the full sheet gravity (dashed curve), 
the sheet monopole approximation (dotted curve), and the softened monopole treatment using /o 
= 1/2 (solid curve). For a fair comparison, all three force profiles are calculated using the same 
surface density and velocity as obtained in §5.1. 
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Fig. 3. — Reduced fluid fields for the softened monopole solutions. The velocity field is shown in the 
lower panel and the density field is shown in the upper panel. In each case, the curves correspond to 
varying degrees of softening, as determined by the zeroeth order flux to mass ratio: /q = 1 (solid), 
/o = 0.75 (long dashes), /o = 0.50 (dashes), /o = 0.25 (dots), and /o = 0.0 (dot-dashed). 
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of the two types of models. Regarding this issue, one should remember that the first three entries 
in Table 1 refer to the same sheet model; for this first model type, the full sheet gravity case is 
nominally the most complete (excluding the central regions of the core). Models of the second 
type are characterized by their /o values, which vary from /o = for the unmagnetized spherical 
limit to /o = 1 for the highly flattened, critically magnetized, case. For /q = 1, the disk thickness 
and the "softening" due to it disappear, and the softened monopole becomes identical to the sheet 
monopole, i.e., they both lack the full complement of higher multipoles that is present implicitly 
in the full sheet gravity model. Nevertheless, when /o ^ 1, Figure [2] indicates that the softened 
monopole model mimics well the best characteristics of the full sheet gravity model while adopting 
none of its pathologies near the origin. For astrophysical applications, we thus regard the results 
presented in the lower portion of Table 1 to be more reliable than those of the upper portion. For 
any of these cases, the associated surface density and velocity field may be reasonably approximated 
with the fitting formulae (j68p and ()69p for o"yi(^) and VAiS,)- 



5.3. Attachment of Condensing Core to Common Envelope 



We now address the solution for flux-to-mass ratio /(^, t). Note that it is possible to discuss 
the effects of ambipolar diffusion only for the softened monopoles of §5.2 where we have made sure 
that -F(^) goes linearly to zero as ^ ^ 0. The magnetic tension force oc — /o-F'(0 ^-cts on the ions 
(but not the neutrals) and drives ambipolar diffusion via the term M in equation (j38p . This term 
would be badly divergent at the origin if -F(^) went to a constant there, rather than vanishing as a 
linear function of i^. As it is, however, the integral of is still logarithmically divergent, as already 
discussed in §3.4. To make M even better behaved, without a lot more analysis, which would be 
an onerous investment when weighed against the limited enlightenment such an effort would yield, 
we adopt the simple procedure of modifying M in the following manner. 



Consider anew the derivation of Appendix B but include from the start the entire current: 

dB 



dt 



1 5 / ^ 

-\ — [wBzu) 

W OW 



where 



^ 1 d 

w ow 



wBi 



dB^ dB, 



(80) 



(81) 



AiT'jppi \ dz dvo 

We replace dB^/dz by {B:^ / Zo)sedi^ {z / zq) and dB^/dw by 2TrG^/'^d{f'S)/diu. Thus, where we 
see B^ we need to add —zofodTi/dzu times some coefficient that represents a thickness correction 
factor. In net, the modified form for the diffusion source term can now be written as 



Ml 



1 



d 



MF + S 



1 da 



(82) 



where M and S are correction factors, respectively, to relate B^ to —F and zodB^/dtu to [(1 — 
/q)/(1 + /g )]cj^^(i(T/(i^ approximately at the origin. Although a proper treatment would require a 
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procedure of singular perturbation theory of the type described in §3.4, we bypass such an involved 
treatment by the simple act of choosing M and S to be the same thickness modification factor that 
we used to regularize F, 

M = , ^ = S. (83) 

With this procedure, A/l goes to a constant at ^ = rather than diverging as 1/.^ near the 
origin. With A/l replacing M , we evaluate the integral 

iVo(0= t^LiOdi. (84) 



^0 

The resulting functions for iVo(^) are shown in Figured] for the zeroeth order flux ratios /o = 0.25, 
0.50, and 0.75. We note that the numerical values for A^o(l) ^-I'e 18.6, 31.1, and 77.4 for the cases 
/o = 0.25, 0.50, 0.75, respectively. Furthermore, the values of A''o(oo) are nearly equal to those of 
A^o(l)) with differences of only ~ 2%. We denote the value of A'^o(l) as a function of /o by the 
symbol /(/o) and we rewrite equation (06]) to obtain the following explicit relationship between e 
and /o: 



(1-/0)71 + 2:^ 







(85) 



Note that the correlation of e with /o is extremely sensitive. Specifically, in order to obtain flux 
to mass ratios /o = 0.75, 0.50, and 0.25, the required values of e are 0.0112, 0.157, and 2.73, 
respectively. We cannot obtain a consistent approximation with e ^ 1 in this formulation when 
/o becomes too low. Conversely, we may say that for standard values of e (~ 0.18), runaway core 
condensation occurs typically when /o reaches 0.5, as indicated by simulations performed by many 
different groups and under very different assumptions (e.g., Nakano 1979, Lizano & Shu 1989, Basu 
& Mouschovias 1994). Values of /o appreciably smaller than 0.5 requires anomalously large values 
of e, under which the condensation problem becomes highly dynamic and can hardly be considered 
diffusive. Such cases, if they exist, need alternative treatments to the one given here. As indicated 
by equation ([85]) . the limiting cases /o = (spherical unmagnetized core) and /o = 1 (completely 
flattened, critically magnetized core) are special in that the first cannot sensibly attach onto a 
critically magnetized envelope unless e is infinite, whereas the second loses flux from its central 
regions and becomes inconsistent with the assumption that /o = 1, the same value as the common 
envelope, unless e is zero. 

Although our derivation of the important equation (|85p was carried out for a very specific 
model, we believe that the result is robust. Indeed, the result is almost given by dimensional 
analysis, except we shall do it through a dimensionless argument. Consider the vector induction 
equation with ambipolar diffusion as it is given by equation (27.12) of Shu (1992), 

- + V X (B X u) = V X I X [B X (V X B)|| . (86) 
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where we have speciaHzed to the ionization law pi = CiocaiP 



(Appendix B). When apphed to 



the core-formation problem, the right-hand side (RHS) of equation (|86p is proportional to three 
powers of B and inversely to the product 7Ciocai- In dimensionless form, the RHS is proportional 
to e/g in the central regions of a condensing core. The left-hand side (LHS) measures the distance 
(in time or distance divided by velocity) that the magnetic field in the core has to travel to get 
from some starting envelope or boundary value to the central value. In proper dimensionless form, 
this distance is 1 — /q. This task is accomplished at the rate oc e/q on the RHS; therefore, when 
the LHS equals the RHS, we have (1 — /o) oc e/g. An order unity quantity on the LHS cannot 
equal an order e quantity on the RHS unless there is a relatively large proportionality factor on 
the RHS. Although this factor might depend on /o (because of the specifics of the model), the 
dependence should be fairly slow since the major dependences should be captured by our scaling 
arguments. The proportionality factor /(/o)/y^1 -|- 2/q in equation ([85|) has precisely these two 
qualities: (1) it is relatively large, and (2) it is relatively constant as a function of /q. Thus, 
apart from minor quibbles about exactly what function /(/o)/\/l + 2/q should be, the relationship 
derived as equation (jSSh is insensitive to the details of geometry, or whether ambipolar diffusion is 
primarily driven by pressure gradients or magnetic tension, etc. 

The physics behind why it is difficult to drive the central flux-to-mass ratio /o to low values is 
now obvious. As ambipolar diffusion occurs and |B| decreases, the rate of diffusion, proportional 
not only to e but also to |Bp, slows down appreciably. It thus becomes increasingly difficult to 
make |B|, relative to S, which is the proper comparison field for |B|, even smaller. As a result, 
many condensing cores get stuck around /o ~ 0.5 before gravitational instability takes over and 
the ratio |B|/S oc / oc is swept into regions close to the origin for further adventures in the 
exciting process called star formation (Mouschovias 1976, Shu et al. 2007). 

To see graphically the formal flux-to-mass profiles implied by our models, we rewrite equation 
(j4ip in the present notation as 



Before runaway condensation occurs, there is not much fluid motion (if we ignore the presence of 
turbulence), and every region can be connected to a nearly static common envelope by character- 
istics from the past, so /(l,ri) = 1. For ^ > 1, Nq{(,) is essentially equal to A''o(l) and we have 
/(^,r) = 1, which identifies such regions as corresponding to the common envelope. For ^ < 1, 
^o{C) < A^o(l) according to Figured! and /(^, r) has a value less than 1, i.e., the common envelope 
makes a transition to a core region undergoing ambipolar diffusion with magnetic flux leaking from 
the ^ < 1 core into the common envelope ^ > 1. In particular, for ^ = 0, Nq(^) = 0, so the 
central flux-to-mass value of the core is /(0,r) = /q. The top panel of Figure [5] shows the flux 
profile /(^) = 1 — [(1 — /o)/A^o(l)] [-^0(1) — -^0(0] before the onset of runaway condensation for the 
realistic cases of /o = 0.75 and /o = 0.50. Because the function Nq{^) is computed assuming values 
of cj(^) and w(^) applicable during runaway condensation (see Fig. this profile should really be 
considered appropriate only for the specific time r = 1 demarcating the transition from quasi-static 



/(C,r) = /(l,ri) 



(1 - /o) 



[No{l)-Nom. 



(87) 
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Fig. 4. — The integral No{^) as a function of ^. The three curves correspond to different values of 
the zeroeth order flux to mass ratio: /o = 0.75 (dashes), /o = 0.50 (solid), and /o = 0.25 (dots). 
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Fig. 5. — Flux to mass ratios. Top panel shows the profiles /(^) for /q = 0.50 (solid curve) and /q 
= 0.75 (dashed curve). Bottom panel shows the flux to mass ratio f{m) as a function of enclosed 
mass m(^) for the same cases, i.e., /o = 0.50 (solid curve) and /o = 0.75 (dashes). 
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evolution by ambipolar diffusion to runaway core condensation. 

After runaway condensation occurs, we can no longer set the advective contribution /(l,ri) 
equal to 1 because ^ > 1 at r ^ 1 lies in the future from = 1 at ti ~ 1 where ambipolar diffusion 
has already occurred to modify / from the value of unity. Instead, /(l,ri) becomes the dominant 
variable term. In comparison, the ongoing diffusion term proportional to e becomes increasingly 
negligible as the difference A'^o(l) — Nq{(^) vanishes because all relevant values of w correspond to 
^ = vu/e^/'^atoT > 1. The physical meaning of this result is that approximate field freezing applies 
during the phases of runaway condensation followed by true gravitational collapse, and the flux-to- 
mass ratio plotted as a function of interior mass becomes fixed in the subsequent evolution. Only 
when very high densities are reached is the assumption of field freezing again violated, perhaps 
involving the formation of a circumstellar disk if we had included the effects of core rotation (see 
the discussion in §7). 

The bottom panel of Figure [5] plots the flux profile / versus the reduced and scaled enclosed 
mass m(^) for the same cases depicted in top panel. This plot shows that the central value /o is a 
substantial under-estimate of the average flux-to-mass ratio of the entire core, i.e., most of the core 
has a flux-to-mass ratio closer to unity. Correction for the under-estimate should act to reduce the 
head-start velocities of real molecular cloud cores in comparison with the values given in Table 1. 
In particular, we can expect extended contraction of the magnitude —Voo only where and when / 
is still climbing to unity. Application of this rule to Figure [5] indicates that extended contraction 
might be observed, perhaps, to about 1/3 of the distance to the outer core boundary (about 1/3 
of the enclosed core mass). As mentioned earlier, these estimates can be made more rigorous using 
singular perturbation theory with multiple length scales, where the region discussed in §5 is treated 
by "intermediate asymptotics" with a scale between the large ones of the common envelope and 
the small ones of the runaway condensation that produces a gravomagneto catastrophe as t ^ 0~. 
Although we have not performed such an improved analysis, we hope that the naive treatment of 
this paper elucidates the physical basis of the phenomenon of gravomagneto catastrophe. 

In any case, at the moment of gravomagneto catastrophe, the enclosed mass in physical units 
inside w = Wce is given by 



where the second equality uses the definition of ©o and defines a benchmark mass scale Mbench = 
a'^Wce/G. For typical values of a = 0.20 km/s and zucc = 0.2 pc, for example, Mbench ^ 1.9Mq. 
The range of /o is limited because values of /o < 0.3 imply values of e greater than unity. Over the 
range from 0.3 to (say) 0.9, (roughly, 1 > e > 0.002), the core masses implied by equation ([88|) vary 
from 2.5Mbench to llMbcnch- Note that this range in mass scale, a factor of 4.4, is smaller than the 
observed range of stellar masses. However, the values of and vuce that specify the mass scale 
Mbench can also vary, and the distribution of these parameter values will add additional width to 
the resulting distribution of core masses. Moreover, final stellar masses can be appreciably smaller 
than the core masses at the beginning of dynamical collapse because of various inefficiencies in 
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the actual star-formation process (such as binary formation and magnetocentrifugally driven winds 
leading to bipolar outflows). These variations will also add width to the distribution of stellar 
masses (see Adams &: Fatuzzo 1996 for greater mathematical detail). Taken as a whole, a strength 
of ambipolar diffusion as a core-formation mechanism is that, given plausible variations of a} and 
tUce, it is capable of producing a core-mass distribution wide enough, when the pivotal state is 
reached, to span the likely pre-collapse states for making brown dwarfs to high-mass stars. 



6. SHEET COLLAPSE SOLUTION 

The analysis presented thus far accounts for the formation of centrally condensed molecular 
cloud cores through the process of ambipolar diffusion. This phase corresponds to negative times 
(t < 0) and can be smoothly matched onto collapse solutions at positive times (t > 0). Although 
the collapse portion of the problem has been considered previously (e.g., Shu 1977, Hunter 1977, 
Gain & Shu 1993, Li & Shu 1997, Krasnopolsky & Konigl 2002, Shu et al. 2004, Fatuzzo et al. 
2004) , here we present a brief re-examination of the problem, and find the particular solution that 
matches onto the solution to the pre-catastrophe problem found in §4 and §5. 

The equations of motion for collapse are the continuity equation ([1]) and the force equation 
([2|). For simplicity, in this treatment we assume flux freezing during dynamical collapse until very 
small scales are reached (see Galli & Shu 1993, Galli et al. 2006, and Shu et al. 2006) so that 
g + £ = (1 — /o )g. Also to keep the discussion uncomplicated, and because there is now a true 
physical monopole (the protostar) to keep the inflowing material spatially flat, we consider the sheet 
monopole limit for the gravitational force (see §4). Here the relevant similarity transformation has 
the form 

x = — , S(t37, t) = — — -(t(x) , and u{w ^t^ = av{x) . (89) 

at ZTT&t 

After some algebra, the dimensionless self-similar form of the equations of the motion become 

dv , . da (x — v) , , 

and 

, dv Gn da , „9, iv — x\ , . 

where Gq and /o have the same meaning as before and are taken to be constants. 
Next we apply the adopted scaling transformation 

i = -=, v = -=, and a = a{x) — . (92) 
Vfo vfo v^'o 

The equations of motion then take the forms 



(93) 
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p|=.(i^|.- «-„)], (94) 

where the discriminant is now given by 

V = {^-vf-l. (95) 
In the outer hmit ^ — > oo, the equations of motion allow the asymptotic forms 

(7 = — and ^ = ^00 + ^ — — ■ (96) 

In order for the collapse solution to match onto the prc-catastrophe solution for t < 0, the constants 
A and Voo must be the same of those of §4.3. Notice also that the sign of the correction term in 
the velocity field is different for the two cases, as it should be: The limit ^ — > oo corresponds to 
t ^ from either side of zero, where the solutions must match, and where both solutions have 
velocity v^o (which is negative, since the core is contracting). For large but finite ^, and negative 
times, the condensation solution has a positive correction to the velocity, so that the velocity is 
smaller in magnitude, i.e., it hasn't reached its full head-start speed that it will at the moment of 
gravomagneto catastrophe t = 0. For large but finite ^, and positive times, the collapse solution 
has a negative correction to the velocity, indicating that the fluid is speeding up as it collapses. 

In the inner limit ^ ^ 0, the scaled equations of motion imply that the solutions have the form 
of a free-fall collapse flow, i.e., 

cr = I 1 fl.nrl n = — I 

where 



mn = lim m(0 = constant . (98) 

Finding the constant tuq is the most important result of the numerical procedure, since the 
rest of the solution is then specified by the similarity transformation. In particular, the dimensional 
mass-infall rate M is given by 

M = -^0 = ^^-0. (99) 
Relative to the standard formula, the scaling of the final expression has the following mnemonic: (1) 

1/2 

the relevant velocity to be cubed is the magnetosonic speed Oq a, and (2) the relevant gravitational 
constant is the magnetically diluted value (1 — /o)G. 

The solution for mg is specified by the pair of constants (A, Voo) that determine how the collapse 
solution for t > matches onto the condensation solution for t < 0. Furthermore, all viable pairs 
of boundary values (^4,^00) correspond to states that are overdense {A > 1) and/or with finite 
head-start velocity — Uoo > 0. These solutions thus correspond to the "outer" solutions in the 
nomenclature of Shu (1977) or the generalization to include nonzero starting velocities (Fatuzzo et 



- 32 - 



al. 2004). In any case, for these "outer" solutions the flow does not go through a critical point. 
As a result, one can directly integrate the equations of motion from asymptotically large ^ (where 
the solution matches onto those of the previous section) down to small ^ <C 1 to determine the 
constant rriQ. 

For the sheet monopole solution to the t < evolution, the boundary values are (^, foo) = 
(1.401, —0.4952). Using these starting conditions, the resulting collapse solution for t > is shown 
as the dashed curves in Figure [6l For this case, the inner constant niQ = 1.670, about 70 percent 
larger than the coefficient found by Shu (1977) for the collapse of the critically stable, singular 
isothermal sphere, mo = 0.975. For the softened monopole solution to the t < evolution, the 
boundary values are {A,Voq) = (1.83,-0.684). For this case, the mass-infall constant for t > is 
niQ = 2.85, and the resulting solution is shown as the solid curves in Figured 

We note that the collapse solution presented here is somewhat idealized, even within the class 
of possible self-similar solutions. The collapse flows shown in Figure [6] are calculated from the sheet 
monopole and softened monopole approximations for the gravitational field. For cases that include 
a full calculation of the perturbational gravity and no initial inward velocities (e.g., Li & Shu 1997, 
Krasnopolsky & Konigl 2002), a shock front develops just outside the infall region. Except for 
the region near the shock, which includes the transition between the inner collapsing flow and the 
outer quasi-static region, the solutions with and without shocks are qualitatively and quantitatively 
similar. Specifically, the collapse solutions with monopole gravity (and no shock front) result in a 
reduced point mass mo « 1.3, whereas the case of full gravity solutions (with a shock front) result 
in mo 1.05 (Li & Shu 1997). In the case considered here, however, the t = configurations (at 
the end of the ambipolar condensation phase and the start of the collapse phase) have non-zero 
inward velocities which act to eliminate the critical points in the fiow (e.g., Fatuzzo et al. 2004), 
so that we do not expect shocks near the head of the expansion wave to play a significant role in 
the collapse. 

This treatment also neglects the effects of rotation on collapse. The solutions found here thus 
represent the outer portion of the collapse fiow, and must be matched onto inner solutions that 
include rotation (Cassen & Moosman 1981, Terebey et al. 1984). When the inner portion of the 
outer region approaches ballistic (pressure-free) conditions, this matching can be done seamlessly 
(Shu 1977; Li & Shu 1997; Fatuzzo et al. 2004). For collapse flows that include magnetic fields, 
however, the roles of magnetic braking and magnetorotational instability (MRI) can be important 
(see Allen et al. 2003; Galli et al. 2006; Shu et al. 2006, 2007). The calculations of this paper show 
that if field freezing strictly holds for the collapse phase t > 0, then the value of Aq brought into 
the star plus disk would be typically ~2, which could prevent disk formation by magnetic braking; 
this result has also been found in numerical simulations (e.g., Fromang, Hennebelle, & Teyssier 
2006, Price &: Bate 2007). How circumstellar disks form and evolve thus remains an open question, 
although it appears likely that global MRI in a context of nonzero net flux will play a major role 
(Shu et al. 2007). 
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Fig. 6. — Reduced fields for collapse solution in the sheet monopole (dashed curves) and softened 
monopole (solid curves) approximations. Upper panel shows the reduced and scaled density field 
cr(^). The lower panel shows the solution for the reduced and scaled velocity field v{^). The initial 
conditions for collapse are taken to be those predicted from the condensation calculation (see Fig. 
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7. ILLUSTRATIVE NUMERICAL EXAMPLE 

To give astronomical context to the semi- analytic results of this paper, we next plot the 
evolution given by the softened monople condensation and collapse solutions in dimensional form 
for the case /o = 0.5 and a = 0.2 km s"^ In Fi gures [7] and [8] we show the equatorial volume 
density p{'Dj,0,t) = T,/2zq, plotted here as the number density n = p/2.3m}i, and the equatorial 
inflow velocity —u{'uj, t) as functions of w and t. Figure [7] shows the time evolution for negative 
times t < (starting from t = —1.0 Myr), whereas Figure [8] shows the time evolution for positive 
times t > (out to t = +1.0 Myr). Notice that the evolution of the runaway condensation phase 
for t < 0, as shown by Figure [TJ compares well with the calculations of ambipolar diffusion carried 
out numerically by Basu & Mouschovias (1994). 

The important thing to carry away from Figure [7] is that condensing cores are not observable 
in dense-gas tracers such as NH3, which requires n > 3 x 10^ cm"'^ for excitation, until the cores 
are within several hundred thousand years of gravomagneto catastrophe. If "cores" are defined as 
such by whether they are observable in dense-gas tracers, then their "lifetimes" will be comparable 
to the lifetimes of embedded protostars, also measured in the several hundreds of thousands of 
years. This numerical coincidence results in roughly equal numbers for "starless cores" and "cores 
with embedded stars," with considerable scatter depending on the value of e in the region being 
studied. Similar statistics given by observers, plus the finding that the surface density profiles of 
cores are flat in their central parts, have led to the mistaken criticism that ambipolar diffusion 
seems to work too slowly to account for the observations (e.g., Andre et al. 1996, Ward-Thompson 
et al. 1999). Runaway core condensation, the phase depicted in Figure O does not take long, 
but it is just the last stage of the ambipolar diffusion process. Indeed, it is a stage where not 
much ambipolar diffusion is still going on, with the central flux-to-mass /o being "frozen" in value. 
Prior to this stage, there were slower stages of evolution, lasting maybe an order of magnitude 
longer than 10^ yr (although it is hard to specify when to start the clock for this less definite 
problem), where ambipolar diffusion did work to get the central regions into the runaway state 
(see Fig. [5]). These stages are not well studied by the similarity methods of this paper but have 
been amply treated by many numerical simulations (e.g., Nakano 1979, Lizano & Shu 1989, Basu 
& Mouschovias 1994, Desch & Mouschovias 2001), and they occupy the bulk of the evolution time 
starting from arbitrarily chosen "initial" conditions. 

The fact that the inflow velocities of "extended contraction" are observed to be a significant 
fraction of the sound speed o indicates that once runaway condensation commences, the process is 
fairly rapid. Nevertheless, if the extended inflow velocities typically reach half the magnetosonic 
speed, the core is roughly within 75% of being "magnetohydrostatic." As a result, the oft-repeated 
statement that "star formation is a dynamical process" is an illusion when applied to i < 0, and 
it is a tautology when applied to t > 0. Observations of extended converging flows thus do not 
support the view of star formation as a turbulent dynamic process, nor do they refute the theory of 
ambipolar diffusion as the core formation mechanism. Indeed, as we argue in §8, the magnitude of 
the observed inflow velocity is fully consistent with the predictions of this paper, and thus provides 
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Fig. 7. — Physical solutions for the number density and velocity of a molecular cloud core during its 
condensation (formation) epoch for varying times before the moment of gravomagneto catastrophe 
(which occurs at t = 0). The displayed time levels are spaced logarithmically, e.g., —0.316 Myr 
separates the curves labeled —1.0 Myr and —0.1 Myr in the top panel. The vertical solid lines mark 
the location of the outer core boundary, if it were set by the condition that the density falls to a 
benchmark value n = 1000 cm^'\ where a core would join onto the background molecular cloud. 
The displayed velocity field is not probably trustworthy when one has reached about 1/3 of the 
distance to the outer core boundary. 
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Fig. 8. — Physical solutions for the number density and velocity of a molecular cloud core during 
its collapse epoch for varying times t > after the moment of catastrophe. The five curves in each 
panel correspond to logarithmically spaced time intervals, i.e., the times +0.0316 Myr, +0.1 Myr, 
+0.316 Myr lie between the displayed times of +0.01 Myr and 1 Myr. The vertical solid lines mark 
the location of the outer core boundary, if it were set by the condition that the density falls to a 
benchmark value of n = 1000 cm^'^, where a core would join onto the background molecular cloud. 
The displayed region of extended infall is probably not probably trustworthy when one has reached 
about 1/3 of the distance to the outer core boundary. 
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a strong argument that molecular-cloud core-formation, at least in isolated regions of relatively low 
total mass, is due to ambipolar diffusion. 

For the collapse phase (t > 0), Figure [8] shows that the fluid fields display the usual forms, 
as studied in many previous treatments (e.g, Shu 1977, Hunter 1977, Galli & Shu 1993, Allen 
et al. 2003, Fatuzzo et al. 2004). In particular, as the collapse proceeds, the densities at a 
given radius w become increasingly smaller than in the pivotal state before it. This behavior is 
a consequence of the solution belonging to the family of expansion-wave collapse solutions (Shu 
1977), with the material missing from the inner core collecting in a point-like object (the forming 
star) at the center of the collapse flow. In the new case considered here, the collapse solution at 
positive times matches smoothly onto the condensation solution of negative times. Not only does 
this generalization provide a self-contained picture of core formation and subsequent collapse, but 
it also shows that the starting state for collapse has a nonzero inward velocity, which results in a 
somewhat larger mass infall rate. The trend of decreasing density at a given radius w is thus more 
extensive than in the case of a magnetohydrostatic starting state, because the "piston" is pulled in 
by the head-start velocities more rapidly and in a more continuous manner. 



8. DISCUSSION AND CONCLUDING REMARKS 

In this paper, we have presented a picture of molecular cloud core-formation that envisages 
a process of the slow leakage of magnetic flux from a dense pocket of gas and dust via ambipolar 
diffusion. This process continues until the central regions acquire a mass to flux to ratio Aq = /q^ 
that is too high for continued quasi-hydrostatic support against the self-gravity of the core, and 
the core develops a runaway central density with a power-law profile S cx or p cc at the 
moment of gravomagneto catastrophe t = 0. One of the key findings of this paper is that /o is 
given by the root of equation ([851) : 

= e (100) 

where K^fo) = /(/o)/y^l + 2/q) ~ 25 for values of /o of physical interest, and e is the dimensionless 
rate of ambipolar diffusion given by equation (j20p . This equation shows that it is difficult for 
ambipolar diffusion to produce regions with /q < 0.3 without an anomalously high effective rate 
coefficient e > 1. A more typical outcome is probably /o ~ 0.5, although values of /o close to 
unity might be sustainable in regions with high ionization rates (which lowers e). Such regions will 
be characterized by high surface densities S oc (1 — /o)~^, and therefore large visual extinctions, 
as well as large core masses, and could be one ingredient to the cores that make massive stars. 
Another ingredient could be high gas temperatures or high levels of turbulence. 

During the epoch of core formation, i < 0, extended regions of contraction develop, and extend 
to perhaps one third of the way to the effective boundary where the core joins onto a common 
envelope of the dense clump that surrounds it (see §5.3 and Fig. 5). Significantly, the contraction 
velocities are always beneath the magnetosonic value. A transition to fully dynamic flow is made for 
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t > that corresponds to an inside-out collapse solution, but with a sub-magnetosonic head-start 
velocity and an over-dense outer envelope given to the region by the ambipolar diffusion process in 
the previous epoch t < 0. The main feature of the t > collapse solution is self-similar infall onto a 
growing protostar at a mass infall rate mo[Oy^a]^/ (1 — /q )G, where the dimensionless coefficient mo 
is two or three times larger than the "standard" value of 0.975 (resulting from purely hydrostatic 
states at t = 0). In some sense, this new solution combines the most attractive features of the 
self-similar solutions proposed in previous work (e.g., Larson 1969, Penston 1969, Shu 1977, Hunter 
1977, Fatuzzo et al. 2004) without bearing any of the unncessary baggage. 

The models predict that the head-start velocities are correlated with over-densities. Compared 
to the static singular isothermal sphere (SIS), the sub-magnetosonically inflowing parts of the core 
characterized by volume densities have total over- density factors of A{1 + 2f^)/{l - f^). For 
/o = 0.5 in the softened monopole model where A = 1.83 (see Table 1), this amounts to an overall 
factor of 3. From a near-infrared extinction study of the bipolar outflow source B335, Harvey et 
al. (2001) found the outer portions of its associated core to have a power-law behavior for 
the inferred volume density, with an over-density relative to the SIS of 3 to 5. They interpreted 
their data in terms of an unstable Bonnor-Ebert sphere. Later, Harvey et al. (2003a) found the 
inner portions of B335 to have a density profile consistent with a r~^/^ law, i.e., consistent with the 
detailed modeling of this source as a classic example of inside-out collapse (Zhou et al. 1993, Choi 
et al. 1995, Evans et al. 2005). In a study of the starless globule L694-2, which has strong evidence 
for inward motions, Harvey et al. (2003b) find the outer portion of the core to have a r^^'^ volume 
density profile, steeper than our model predictions, but with an over-density factor in its central 
regions relative to Bonnor-Ebert extrapolations of about 4. While these authors speculate that the 
effect might arise from the core being an prolate object viewed along its long axis, a more satisfying 
and unified interpretation is that B335 is a t > post-catastrophe core with an embedded protostar 
(and bipolar outflow), magnetized at /o ~ 0.5; and that L694-2 is a t < pre-catastrophe starless 
core, magnetized also at /o ~ 0.5. Another indicator of the correctness of this identification is 
that the predicted mass-infall rate M = uioIQq^'^ a]^ / {1 ~ fo)G is roughly consistent with measured 
values in Class sources (for examples with estimates at the extremes, see Ohashi et al. 1997 and 
Furuyu et al. 2006), despite earlier claims that Class protostars would have much higher infall 
rates (e.g., Henriksen et al. 1997). 

The (near) self-similarity of the problem is a particularly attractive feature of the process. 
Given that the central portions of the core are nearly isopedic, i.e., that A = 1// is nearly a spatial 
constant, self-similarity of the collapse solution for t > (but not too much greater!), or even for 
the runaway core-condensation phase for t < (but not too much less!) is perhaps not a surprising 
outcome of nature's tendency to produce power-laws when solutions have to span large dynamic 
ranges in space and time. 

However, the reader could rightfully question whether the /-profiles obtained in Figure [5] are 
not special to the application of self-similarity to a problem - the initial stages of the condensation 
of a cloud core by ambipolar diffusion - that has no good reason, beyond mathematical convenience. 
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to be self-similar. After all, the contraction of typical molecular densities in clumps of ~ 10^ cm~^ 
to early-stage cloud cores with number densities ~ 3 x 10^ cm~^ can hardly be characterized as 
spanning a very large range. One might think that the resulting /-profiles would show considerable 
variation depending on the exact initial state being assumed and boundary conditions being applied. 
And so it must be with very detailed descriptions of such early stages of cloud evolution. But if one 
is pressed for more global trends, there are only so many ways that a function / can monotonically 
go from unity at some large radius to some other value /o, typically 1/2, at some small radius. 
And a self-similar approach to getting such a profile is probably not any worse than some other 
ad hoc prescription. The important features of the picture are not the details of the /-profiles, 
but the global view provided by estimates of the relevant time scales, relationships between e and 
/o, generic stages of the evolution, and the final asymptotic convergence to self-similarity as the 
moment of gravomagneto catastrophe is approached and passed. 

In particular, the solutions presented in this paper are compatible with full ambipolar-diffusion 
calculations that start with marginally subcritical configurations which develop nearly isopedic cen- 
tral cores (with A ^ constant) before proceeding on a path of extended gravitational condensation 
that leads to gravomagneto catastrophe (Nakano 1979, Lizano & Shu 1989, Basu & Mouschovias 
1994). It is particularly significant that both the predicted and observed contraction velocities are 
sub-magnetosonic and arise from the modest over-densities that are left behind in the contracting 
cores as their magnetic support leaks to the common envelope. Thus, the observed head-start 
velocities are an indicator that some slow process like ambipolar diffusion is at work producing 
molecular cloud cores, rather than some more sudden process of the destruction of high levels of 
non-thermal support, such as the dissipation of hypersonic turbulence through shock waves. The 
latter description may still apply, however, in the crowded conditions that characterize high-mass 
star-forming regions. 

The most important lesson of this paper is that the details of the ambipolar diffusion process 
control only the spatial extent of the region of extended, sub-magnetosonic, contraction and the 
timing of the runaway core-condensation that leads to the gravomagneto catastrophe. Provided 
the small parameter e is not strictly zero, gravomagneto catastrophe is the unavoidable fate of a 
lightly- ionized, isolated, molecular-cloud core, as long as Lorentz forces contribute to the support 
against its sclf-gravitation (see the discussion of Lizano & Shu 1989 concerning "failed cores"). 
Given that the inner cores acquire nearly isopedic states with A = 1// ~ constant, the resulting 
density and magnetic field profiles of the resulting runaway condensation steepen into generic power 
laws and are robust. 

This work was initiated during a sabbatical visit of FCA to UCSD; we would like to thank 
the Physics Departments at both the University of Michigan and the University of California, San 
Diego, for making this collaboration possible. Discussions with Mike Cai were very helpful. This 
work was supported through the University of Michigan by the Michigan Center for Theoretical 
Physics; by NASA through the Astrophysics Theory Program (NNG04GK56G0) and the Spitzer 
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A. Magnetic Forces in a Thin Disk 

In this Appendix, spurred by a correction pointed out by Mike Cai (2007, private communi- 
cation), we revisit the derivation given by Shu Sz Li (1997) for the magnetic forces in a thin disk. 
The Lorentz force per unit volume with only axisymmetric poloidal fields is given by 

-1(V X B) X B = - ^) (B^e^ - B^e,) . (Al) 

47r 47r \ ovj oz J 

If we integrate over z in a thin disk of effective thickness 2zq <^ w with = —B^ being, 
respectively, the radial magnetic field at the upper and lower surface of the disk and with B^ being 
continuous across the midplane, the force per unit area in the radial direction is given by 



27r dw V 47r 



B,B+ d (B'lzo 



(A2) 



We recognize the first term as what Shu k, Li refer to as the force per unit area due to magnetic 
tension, but the second term is only the negative radial gradient of the vertically integrated magnetic 

pressure due to B'^/Stt and not {B^ + B'^)/8tt. The reason is that the so-called "magnetic-tension 
term" also contains a small piece of the magnetic pressure, in fact, exactly the integral of B'^/Stt 
over the disk thickness. Nevertheless, for simplicity we shall continue to refer to the two terms as 
"magnetic tension" and "magnetic pressure." 

Vertical hydrostatic equilibrium for a magnetized isothermal gas in its own vertical gravita- 
tional field requires 

where the last term is the dominant term for the magnetic force per unit volume in the z direction 
according to equation (Al). In the above, U is the self-gravitational potential of the gas and satisfies 
the local Poisson's equation: 

-Q^ = ^^Gp. (A4) 
The substitution of the above into equation (A3) allows us to integrate once: 

^ ^f^V + «V + |^ = C(^), (A5) 



SvrG \dz J Stt 

where C{'ud) is a constant for fixed w (and diffusion time). We evaluate C at the upper disk surface 
where dU /dz = 27rGS, /j = 0, and B^ = B^; and we do the same at the disk mid-plane where 
dU/dz = 0, p = Ti/2zQ (being a definition of zq), and = 0. Setting equal the two expressions 
for the "constant" total pressure, we obtain 

2zQ 2 Stt ^ ^ 
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For an isopedic singular isothermal disk, which is what the inner parts of molecular cloud cores 
become at the moment of gravomagneto catastrophe, 

B^ = B, = , (A7) 

with A equal to a constant. Equation (A5) now becomes the second relation of equation ([8]). In 
the same limit, we have 

B^zo ( a^H 



which shows that the sum of gas pressure force in the radial direction with the second term in 
equation (Al) equals —a'^d{@'S)/dw with G given by the first relation of equation ([8]) rather than 
by the expression e = (3 + A2)/(1 + A^) from the analysis of Shu & Li (1997). The difference (2 
versus 3 in the sum of the numerator) arises because the latter authors mistakenly included the 
contribution of (B^)'^ /8tt into the computation of the magnetic "pressure" force, which duplicates 
a small piece already included in the first term of equation (Al). 

For later reference, we consider the vertical structure if we make the assumption that the 
current density is proportional to the volume density, i.e., 

dB^ _ ^2p 

The substitution of equation (A8) into equation (A5) yields a differential equation which we may 
write as 

1 f9UY^^r2B^^Bl^.^^,^iB^^ j^^j 



SvrG \dz J 2Bt dz Svr 2 Stt 
The solutions for equations (A4), (AS), and (A9) read 

B^ = 5+ tanh(z/zo), P= -^sech^{z/zo), = -27rGi; tanh(z/zo), (AlO) 

ZZq OZ 

if we set B^ = B^ = 2t:G^/'^Yi/ \. Equation (AlO) represents the isopedically magnetized version of 
Camm (1950) and Spitzer's (1955) solution for the stratified isothermal disk. We note in passing, 
however, that while equation (A7) is an acceptable approximation in the disk proper, it must fail 
on the z-axis where B^ = from symmetry considerations, but B^ is not only nonzero but achieves 
a maximum value there. 



B. Derivation of the Diffusion Constant 



In this Appendix, we present a brief discussion of the derivation of the diffusion coefficient in 
the magnetic diffusion equation. We start with the z component of the induction equation including 
ambipolar diffusion 



dt 



1 a Id \ wBj dB^ 

wow wow liTT'yppi OZ 



(Bl) 



where we have kept only the largest term on the right hand side for a highly flattened core. Here, 
is the z component of the field and B^ denotes the vj component. After multiplying the equation 
by p and integrating over z, we obtain the form 



dB, I d , 

H ^ i^B^u) 

at w ow 



1 d 


\wBl 


dzdB^- 


1 d 


'wBlBl{2z,r^' 


w dw 




Loo Pi dz 


w dw 


27r7CSi/2 



(B2) 



where we have defined 



1 

C 



SV2 



dz dBr, 



(2zo)V2B+7o ft dz 



sech(z/zo) dz 

Clocal 



(B3) 



if we assume equation (AlO) to hold and express the ion abundance by the local relation pi = 
CiocaiP^''^ (Shu 1992). In the present application, Ciocai is a constant if cosmic-rays provide the 
dominant source of ionization, but it quickly climbs to much larger values near the surfaces of 
molecular clouds because of the ultraviolet ionization of elements like carbon (McKee 1989). For 
Cioeai = const, we have C = (2/7r)Ciocai = 2.0 x 10"^^ cm^^/^ gi/2 ^gj^^ ^992)_ 



C. The Velocity Function is Monotonic 

In this Appendix, we argue that the reduced and scaled velocity field u is a monotonic function 
of ^ for the regime of interest. Since t; = at = (the inner boundary condition), the monotonicity 
of V implies that the solution must have a nonzero velocity at large ^. In physical terms, this finding 
implies that starless cores are predicted to have nonzero velocities, even before the collapse phase 
begins (as observed). Of course, the numerical integration of the equations of motion implies 
nonzero values of Voo- In this Appendix, however, we analytically show that this property must 
always hold. 

In order to prove this assertion, at least in the context of the approximations of this paper, it 
is sufficient to show that the right hand side of equation (j58p is never equal to zero except at the 
critical point (where the discriminant V changes sign). First, we define the ancillary function 

P{i)^<y{i + v). (CI) 

The right hand side of equation (|58p will be zero if and only if P{C) = 1- Further, we know that 
P(^) = 1 at the critical point. Next we show that -P(^) is monotonic. Differentiating P with respect 
to ^, we obtain 

dP ,^ .da / dv\ 

_ = K + + + (C2) 

Using the equations of motion (eqs. [59] and [58]), this result simplifies to the form 
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Since the density a and the coordinate ^ are always positive, and since we are interested in con- 
tracting solutions where v is negative, the right hand side of this equation, and hence dP/dS, is 
positive. As a result, the function -P(^) is a monotonically increasing function of the variable ^. 

Since P{(,*) = 1 and P{£,) is monotonic, it follows that P < 1 for all ^ < and P > 1 for all 
^ > It then follows that the right hand side of equation (1581) is positive for < ^* and negative 
for ^ > 1^*. Since the discriminant has the opposite behavior, T> < for ^ < ^* and T> > for 
^ > it follows that —dv/d^ > for all values of ^. Thus, v is a monotonic function of as 
claimed. 



D. Generalized Approach to Crossing Critical Lines 

In this Appendix, we record the procedure needed to solve the posed problem when the force 
integral p4p has a general form written as 

F.-A({)!fi = -A({)«i^. (Dl) 

The first equality is true by definition, i.e., we define the function A(^) to be the ratio of the true 
force to that given by the monopole approximation. The second equality follows from the continuity 
equation. THe monopole approximation corresponds to the simplest case A(^) = 1. 

With the introduction of the correction function A(^), we can find the values of the fluid flelds 
and their derivatives at the critical points. Speciflcally, for critical point we find 

^(^*) = T7TT and vi^,) = 1 - ■ (D2) 

Using the same expansion around the critical point as before (see eq. j61j). we find the derivatives 
of the fluid flelds at the critical point, i.e., 

v, = -^±±- [(e, - 1)2 + 1 - 2e.A:/A,] , (D3) 

and 

= {e* - 2 + 2e*A'JA, T [(e* - l)' + l - 2^.A'JA,] . (D4) 

In these expressions. A* = A(^*) and A'^ = dA/d^{^^). 



E. Multipole Approach to Effective Sheet Gravity 

In this Appendix, we consider the multipole approach to the evaluation of the force equation 
(see, e.g., Li & Shu 1997), which we write as 

dU 

^(0 = -^, where ^^(0 = Hoi^ , vMri) r,dr, , (El) 
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with Hq being the classical Poisson kernel for a self-gravitating axisymmetric sheet, 

dip 



_|_ ^2 _ 2^j^ COS ip 



(E2) 



We now use the well-known formula for the spectral expansion of the inverse separation distance 
between a field point at ^ and a source point at 77 separated by an angle ip (see, e.g., eq. 3.41 of 
Jackson 1962): 



00 f 

7] 



if + _ 2^r? cos (^)- V2 = J2 7^^!^^ (cos for r? < 



for 



(E3a) 
(E3b) 



where Pi{fJ,) are the Legendre polynomials of order £. 

For an axisymmetric surface density distribution (7{r]), we now have 



2e 



where 



P2i (cos ip)dip, (E4a) 



and 



U<= ri^'aiv)vdv, U>= r?-(2^+i)a(r?) r?dr?. 
Jo if 



(E4b) 



Our sum extends over only even values of £ because the coefficients vanish for odd i. The 
numerical values of C2e are all positive; according to (an equivalent formula by) Gradshteyn & 
Ryzhik (1980, 7.222): 



C2e 



(2£- 1)!! 
(2^)!! 



1 2 



Thus, C2i = 1, 1/4, 9/64, 25/256, 1225/16384, 3669/65536, etc, for 2£ = 0, 2, 4, 6, 8, 10, etc, 
declining slowly only as ~ l/2£ for large £. We refer to C/^ and [/^ as, respectively, the interior (or 
inner) and the exterior (or outer) multipole moment of order 2£ with 2£ = 0, 2, etc, corresponding 
to the monopole, quadrupole, etc. In any case, if we now carry out the differentiation indicated in 
equation (El), we get 



e=o 



(E5) 



where we have used the fact that each multipole order cancels in pairs if we differentiate the moments 
rather than the powers of ^. Note that only the interior monopole moment term oc Uq survives 
this differentiation because 2£C/^ equals zero when £ = 0. This well-known result is fortunate since 
Uq is formally logarithmically divergent if (7(77) — A/ry at large 77. 
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If cj(r?) = A/r] for all 77 from to 00, then U< = ^2^+V(^ + 1) and Uf = i'"^^ /2l (except for 
^ = 0), and all multipoles cancel in pairs in F except for £ = 0. As is well-kown, the radial force 
field for a perfect singular isothermal disk (SID) is given solely by the interior monopole. When 
<T(ry) departs from the ideal SID state, say by becoming a constant ctq in the central regions, then 
U21 ~ '7oC'^^^^/(2^ + 2) for small ^, so each interior multipole contributes a constant term to F{^) 
at small ^. But the exterior multipole moments will now contribute terms that are larger, and 
opposite in sign, to their interior multipole counterparts. Thus, the effect of including multipoles 
reduces the inward force of F[^) at small ^ relative to the monopole contribution for given A. 

Computing F{^) by equation (j74p is equivalent to summing the infinite set of multipole con- 
tributions. In either case, we can compute the correction function A(^) of Appendix D as 

m = (E9) 

m 

where m = Uq. If the correction function A(^) were known, then one could find the solution using 
the same procedure as before (in §4 for the monopole solution): Guess the value of the critical 
point, move inward from the working estimate of ^* using the results of §6.1, and integrate inward 
to the origin. Then adjust the value of the critical point and iterate until the inner boundary 
conditions are satisfied. After finding the critical point, one further integration of the equations of 
motion (both inward to the origin and outward to large ^) then determines the solution. In this 
case, however, we do not know the function A(^) and its form depends on the solution for <t{^) 
that we are trying to find. As a result, we must use another iterative scheme: We first estimate 
(guess) the form of the function A(^), and then calculate the (approximate) solution according to 
the previous procedure. With this approximation to a, we can evaluate the integrals in equations 
(E4a) and (E4b) to find a new estimate for the correction function A(^). We then iterate this 
procedure until the solution is obtained. 



F. Full Effective Gravity of Unfiattened Core 

In this Appendix, we consider the properties of the full gravity of an incompletely flattened 
core. This can be carried out by replacing the kernel 'Ho(C,^) by the weighted-average of the 
product of the volume densities at the source and field points of the 3-D Poisson integral (see eq. 
A9): 

1 f^^ dC dC sech2[C/Co(0]sech2[C/Co(r/)] 

^i?, '/; = -■ 

where 



CoK)- 



2(1 - fl 
1 + 2/0' 



^(0 



(F2) 



It is trivial to show that d7i{0,rj)/d^ = 0. Physically, an axisymmetric magnetized disk with finite 
thickness cannot exert a net radial force at its center if its volume density and current density are 
regular there. 
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Although it is possible to develop a multipole expansion procedure for equation (Fl), the 
resulting analysis would be quite involved. For simplicity, therefore, we are content to adopt the 
alternative treatment of §5.2 designed to give a physical assessment of the influence of finite disk 
thickness in the "softened monopole" approximation. 
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Table 1: Parameters for Diffusion Epoch Solutions 



Model 








A 


Monopolc 


1.294 


5.03 


0.495 


1.40 


Quadrupole 


1.407 


3.12 


0.732 


1.68 


Full Sheet Gravity 


1.714 


1.98 


1.67 


3.43 


/o = 0.0 


1.467 


3.08 


0.833 


2.35 


/o = 0.25 


1.455 


3.04 


0.815 


2.21 


/o = 0.50 


1.393 


3.22 


0.684 


1.83 


/o = 0.75 


1.328 


3.85 


0.556 


1.52 


/o = 1.00 


1.294 


5.03 


0.495 


1.40 



